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In this paper the two dimensional extended Hubbard model with intersite magnetic Ising-like 
interaction in the atomic limit is analyzed by means of the classical Monte Carlo method in the 
grand canonical ensemble. Such an effective simple model could describe behavior of insulating 
(anti)ferromagnets. In the model considered the Coulomb interaction (U) is on-site and the mag¬ 
netic interactions in ^-direction (J > 0, antiferromagnetic) are restricted to nearest-neighbors. Sim¬ 
ulations of the model have been performed on a square lattice consisting of A/" = L x L = 400 sites 
(L = 20) in order to obtain the full phase diagram for U/{4J) = 1. Results obtained for on-site 
repulsion (U > 0) show that, apart from homogeneous non-ordered (NO) and ordered magnetic (an- 
tiferromagnetic, AF) phases, there is also a region of phase separation (PS: AF/NO) occurrence. 
We present a phase diagram as well as some thermodynamic properties of the model for the case of 
f//(4J) = 1 (and arbitrary chemical potential and arbitrary electron concentration). The AF-NO 
transition can be second-order as well as first-order and the tricritical point occurs on the diagram. 
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I. INTRODUCTION 

Since its introduction in 1963 [1] the Hubbard model 
has found applications in many various systems. Despite 
half-century research on this model it still holds a num¬ 
ber of open questions. This report focuses on the atomic 
limit of the extended Hubbard model, in which we restrict 
ourselves to the case of the zero-bandwidth limit (t = 0). 
with added magnetic interactions of the Ising-type be¬ 
tween electrons. Such a simple model can be used for 
describing behavior of insulating magnets. The hamilto- 
nian of the discussed model has the following form: 

H = u'^hi-fhii + 2J s|s| - hi 
i (ij) i 

where U is the on-site density interaction, J is 2 :- 
component of the intersite magnetic exchange interac¬ 
tion, /i is chemical potential, and restricts the sum¬ 

mation to nearest-neighbor sites, rii = fii^ + is total 
electron number on site i, whereas si = (1/2) (n^^^ — n^) 
is 2 :-component of total spin at i site, fiia- = is 

electron number with spin a on site i, where and Qct 
denote the creation and annihilation operators, respec¬ 
tively, of an electron with spin a =t, i at the site i. The 
electron concentration n is defined as n = {1/N) 
where N is the total number of sites. 
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The rigorous ground state results for this model have 
been found in the case of a d = 1 chain 0 E] and 
2 <d < +00 case 0 0. The exact results for finite 
temperature have been also obtained [3] for (i = 1 chain 
(an absence of long-range order at T > 0). Within the 
variational approach the model has been analyzed for 
half-filing (n = 1) [HI E] as well as for arbitrary electron 
concentration 0 < n < 2 [H [9] (these results are rigor¬ 
ous in the limit of infinite dimensions d +oc). Our 
preliminary Monte Carlo (MC) results have been pre¬ 
sented in m for strong on-site repulsion {U / A J = 1,10 
and L = 10). In this paper we investigate in details 
the phase diagram and thermodynamic properties of the 
model for arbitrary electron concentration n < 1 and ar¬ 
bitrary chemical potential ft < 0 {]! = fi — U/2) in the 
whole range of temperatures for a specific repulsive value 
of the on-site interaction parameter U/(AJ) = 1 (and 
L = 20). The corresponding results for n > 1 and /i > 0 
are obvious because of the electron-hole symmetry of the 
model on alternate lattices. 


II. MONTE CARLO SIMULATION DETAILS 

The Monte Carlo simulations for the model described 
above have been done at finite temperatures T > 0 us¬ 
ing grand canonical ensemble on two dimensional {d = 2) 
square (SQ) lattice with number of neighbors z = A. One 
could depict such approach as adsorption of electron gas 
on a lattice. Our simulations use a local update method 
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FIG. 1. Magnetization tuq as a function of temperature 
ksT/AJ for different system sizes L = 4, 6, 10, 20, 30 (as 
labelled) for U/AJ — 1 and Jl/AJ — —0.22. 

m determined by elementary “runs”: (i) a particle can 
transfer to another site, (ii) it can be adsorbed on a lat¬ 
tice or (iii) it can be removed from the lattice. These 
“runs” are usually called move^ create and destroy of par¬ 
ticle procedures in MC simulations. Every Monte Carlo 
step (MCS) consists of each of these “runs” performed 
N = L X L times. In our simulations the number of MCS 
is 10 ^ with a quarter of them being spent on thermaliza- 
tion, which is necessary to avoid results heavily influ¬ 
enced by the starting point of the simulations. There are 
also many cluster update algorithms, but due to chemical 
potential term in the hamiltonian they cannot be imple¬ 
mented here. More details on the simulation method can 
be found in m 

Simulations provide data for temperature and chemical 
potential dependencies of various thermodynamic vari¬ 
ables. The variables of the particular interest are: stag¬ 
gered magnetization tuq = {rriA — tub)I‘2, magnetic sus¬ 
ceptibility XruQ = {{mg) - {mQ)^)/{TN), and specific 

heat c = ((E^) - {E)‘^)/{T‘^N) {E = (H)). Because anti¬ 
ferromagnetic interactions (J > 0 ) are studied, staggered 
magnetization mg is an order parameter in the model 
considered, which is calculated as a difference of magne¬ 
tization of sublattices A and B {rria = {2/N) 
a = A,B). In the antiferromagnetic (AF) phase stag¬ 
gered magnetization is non-zero {mg ^ 0 ), whereas in 
the non-ordered (NO) phase mg = 0 . 

III. RESULTS AND DISCUSSION (U/(4J) = 1) 
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FIG. 2. Electron concentration n (red), magnetization 
m (blue), magnetic susceptibility Xmq (green) and specific 
heat c (violet) as a function of chemical potential y/4J 
for ksT/AJ = 0.1S (a) and /csT/4J = 0.36 (b); and as a 
function of temperature ksT/AJ for /i/4J = — 0.67 (c) and 
^/4J = -0.18 (d) (all for U/AJ = 1 and L = 20). 


The transitions in finite systems are not sharp and 

the finite-size effect on the order parameter mg is ob- T = 20 yieids an essentiai change in the resuits, a further 

served in the resuits of the MC simuiations. In the NO increase of L does not make the transition sharper and 

phase mg is iarger than zero {mg 7 ^ 0) even above the greatiy increases simuiation time. Thus system size of 

AF-NO transition temperature. The temperature de- N = L x L = 400 has been chosen and aii further resuits 

pendence of magnetization mg for different SQ iattice are for L = 20 uniess said otherwise, 
sizes is shown in Fig. Whiie a change from L = 10 to For a given set of modei parameters temperature and 












































3 



FIG. 3. (a) The ksT/AJ vs. /i/4J phase diagram for 

t//4J = 1. T denotes a tricritical point. Dashed-dotted lines 
indicate diagonals labeled as 20 and 30. (b), (c) Staggered 
magnetization mq for heating (red) and cooling (blue) runs 
for diagonal index D = 20 (first-order transition) (b) and 
D = 30 (second-order transition) (c) as a function of tem¬ 
perature ksT/AJ (all for U/AJ — 1 and L — 20). 


chemical potential dependencies of thermodynamic prop¬ 
erties (n, mg, XmQi and c) have been obtained as illus¬ 
trated on Fig. A location of critical points is done 
by analysis of mg, Xtuq and c. The traversal of the 
boundary between two phases (AF and NO phases) is 
usually connected with a substantial change of magne¬ 
tization mg. However, because the finite size effects is 
observed in the dependence of mg, a more precise loca¬ 
tion of the critical point (i.e. AF-NO transition) is de¬ 
termined by the discontinuity of magnetic susceptibility 
XrriQ as well as a peak in c. 

Instead of running simulations with fixed ksT/AJ or 


the simulations with both these thermodynamic 
parameters changing have been also performed. In such 
an approach the system moves diagonally in (/i, T) plane, 
approaching the phase transition area almost perpen¬ 
dicular, which allows for more precise determination of 
the AF-NO transition temperature. We introduce def¬ 
inition of a diagonal {D)^ where a diagonal labelled as 
0 is perpendicular to /c^T/dJ-axis and a diagonal la¬ 
beled as 100 is parallel to it. For each step on the di¬ 
agonal D the parameters change as ATjj = I) AT/100, 
ApD = (100 — T)A/i/100, where D is diagonal number, 
and AT, AjU correspond to temperature (ksTlAJ) and 
chemical potential (ft/AJ) fixed steps on the respective 
axes. The final ksT/AJ vs. Jl/AJ phase diagram for 
UjAJ = 1 is shown in Fig.j^a). In Fig.j^a) there are also 
indicated the diagonals with indexes T = 20 and D = 30. 
Staggered magnetization mg for heating with removal of 
particles and cooling with addition of particles for diag¬ 
onal index 20 and 30 as a function of temperature for 
UfAJ = 1 are shown in Figs.j^b) and (c), respectively. 

The AF-NO transition temperatures increase mono- 
tonically with decreasing |/i|. The maximum of the 
transition temperature is located at p = 0 (n = 1). For 
U/AJ = 1 and n = 1 the critical temperature is equal to 
kBT/AJ-0.A7 (T = 20). 

To determine the location of the tricritical point T 
two simulation runs have been done for each diagonal, 
one starting at (0,0) in the p/AJ-ksT/AJ plane and 
the other starting at the maximum point for the given 
diagonal and descending to (0,0). This corresponds to 
heating and cooling processes, respectively. A position 
of the T point can be estimated by comparing magne¬ 
tization data at the point of the phase transition for 
those two simulations. For second-order phase transition 
both magnetization curves should be almost identical 
(Fig. ^c)), while first-order phase transition is character¬ 
ized by the hysteresis (Fig. |^b)). Thus, a point, where 
the hysteresis is collapsing into single curve, is a tricriti¬ 
cal point. For U/AJ = 1 (L = 20) the T point is located 
at ksT/AJ = 0.205 ± 0.005 and p/AJ = -0.895 ± 0.007 
(n = 0.52 ±0.2). 

With simulations done for fixed p, it is possible to 
obtain phase diagrams as a function of n by determin¬ 
ing electron density above (n_) and below (n+) the 
AF-NO phase transition (for fixed p). The first-order 
AF-NO boundary for fixed p splits into two boundaries 
(i.e. PS-AF and PS-NO) for fixed n. For U/AJ = 1 the 
ksT/AJ vs. n phase diagram is presented in Fig. For 
temperatures above the tricritical point T, the AF and 
NO phases are separated by second-order line. At lower 
temperatures, below this point, there is a phase sepa¬ 
rated (PS: AF/NO) state. The PS state is a coexistence 
of two (AF and NO) homogeneous phases. 

Comparing the results presented in this paper with pre¬ 
liminary MC simulations from previous work m of our 
group, the general improvement of their quality is clearly 
seen. Increasing the SQ lattice size as well as simulat¬ 
ing along diagonals gives more accurate values of critical 
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FIG. 4. The ksT/AJ vs. n phase diagram for Uj^J = 1 (and 
L = 20). T indicates a tricritical point. 


temperatures. Moreover, performing heating and cooling 
simulation runs yields much more precise location of the 
tricritical point. 


IV. FINAL COMMENTS 

The results presented in this paper are in a good qual¬ 
itative agreement with those obtained within the varia¬ 
tional approach (VA) involving the mean-field approxi¬ 
mation for intersite interactions dun], which is exact in 
d +OC (L oc, z ^ oo). Our MC simulations have 
been performed for d = 2 SQ lattice. In this case the 
VA is much less reliable. It largely overestimates critical 
temperatures, e.g. ksT/AJ 0.68 at n = 1, and yields 
a quite different location of the T point: ksT/AJ = 1/3, 
Jl/AJ —0.96, n 'Zi 0.34 (for U/AJ = 1) O [10]. While 
finite-size effects do not pose a big problem, long thermal- 
ization time at low temperatures prevents from obtaining 
results near the ground state. At such low temperatures 
the probability of an electron to change its state is min¬ 
imal, so the system has little chance of escaping false 


energy minima. A solution of this problem is running 
simulations for a very long time, or collecting results from 
different “runs” for various starting states. 

It is important to mention that in the absence of an ex¬ 
ternal magnetic field the antiferromagnetic (J > 0) inter¬ 
actions are simply mapped onto the ferromagnetic ones 
(J < 0) by redefining the spin direction in one sublattice 
in alternate lattices decomposed into two interpenetrat¬ 
ing sublattices. Thus, our results obtained in this paper 
are still valid for J < 0 if tuq m = {rriA + mB)l‘2 and 

The analysis of effects of finite band-width {t ^ 0) is 
a very important problem. However, because of the com¬ 
plexity of such model only few results are known beyond 
weak coupling regime or away half-filling (ME]- For in¬ 
stance, the presence of the hopping term ^ Ujcl^Cja- 
breaks a symmetry between J < 0 and J > 0 cases nu¬ 
ns]. The detailed analysis and discussion on this topic is 
left for future study. 

The competition between magnetism and supercon¬ 
ductivity [2QH26] in atomic limit of the extended Hubbard 
models is a very interesting topic. Some results concern¬ 
ing the interplay of magnetic interactions with the pair 
hopping term have been presented in EOj. Moreover the 
interplay between various magnetic and charge orderings 
[?7f[33] has been also analysed [9l|34|. 
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